Stochastic analysis of surface roughness 
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Abstract 

For the characterization of surface height profiles we present a new 
stochastic approach which is based on the theory of Markov processes. 
With this analysis we achieve a characterization of the complexity of the 
surface roughness by means of a Fokker-Planck or Langevin equation, 
providing the complete stochastic information of multiscale joint proba- 
bilities. The method was applied to different road surface profiles which 
were measured with high resolution. Evidence of Markov properties is 
shown. Estimations for the parameters of the Fokker-Planck equation are 
based on pure, parameter free data analysis. 

1 Introduction 

The complexity of rough surfaces is subject of a large variety of investigations in 
different fields of science JIBGIll . Physical and chemical properties of surfaces 
and interfaces are to a significant degree determined by their topographic struc- 
ture. As one example, the influence of surface roughness on boundary layer flows 
is discussed in turbulence research (cf. 5,6 ) and in the atmospheric (cf. J7HB1) 
and oceanographic sciences (cf. [5]). A comprehensive characterization of the 
topography is of vital interest for deposited, polished or otherwise processed sur- 
faces. Therefore, to give a second example, current roughness analysis methods 
find practical application for the characterization of polished surfaces [10111112] . 

In the context of industrial and engineerig applications common roughness 
measures are based on standardized procedures as the maximum height differ- 
ence R z and the mean value of absolute heights R a These rather simple 
measures clearly cannot completely characterize the complexity of roughness. 
This is also confirmed by the existence of a large amount of additional measures, 
each describing a very special feature of a surface. 

In the physical and other sciences a common local measure of roughness is 
the rms surface width w r (x) — ((h(x) — h) 2 )V 2 , where h{x) is the measured 

* E-mail: matthias . waechterOuni-oldenburg . de 
^E-mail: peinke@uni-oldenburg.de 



height at point x, ( • ) r denotes the average over an interval of length r around 
the selected point x, and h the mean value of h(x) in that interval. Advantages of 
w r (x) are a scale dependent definition as well as clear physical and stochastical 
meanings. Among the techniques used to characterize scale dependent surface 
roughness probably the most prominent ones are the concepts of self-affinity 
and mult i- affinity, where the multi-affine /(a) spectrum has been regarded as 
the most complete characterization of a surface 14, 15, 2, 3 . From a stochastic 
point of view we want to point out two important properties of multiaffinity: 
(a) the ensemble average (w r ) must obey a scaling law (w") ~ r^ a , and (b) 
the statistics of w r {x) are investigated on distinct length scales, thus possible 
correlations between w r (x) and w r >(x) on different scales r, r' are not examined. 

The method we are proposing in this contribution is based on stochastic 
processes which should grasp the scale dependency of surface roughness in a 
most general way. No scaling feature is explicitly required, and especially the 
correlations between different scales r and r' are investigated. To this end we 
present a systematic procedure how the explicit form of a stochastic process 
for the r-evolution of a roughness measure similar to w r (x) can be extracted 
directly from the measured surface topography. This stochastic approach has 
turned out to be a promising tool also for other systems with scale dependent 
complexity like turbulence |16I17| and financial data ^8] ■ I n this letter we focus 
on the complexity of rough surfaces. As a specialized example we have picked 
out the applied problem of characterizing road surfaces. These are an essential 
component of current transportation and thus represent a class of non-idealized 
and widely used surfaces. We claim that the applicability of our method to this 
class of surfaces indicates its general relevance for improved surface characteri- 
zation. 

A collection of road surface data was measured which served previously as 
an empirical basis for the prediction of vibrational stress on bicycle riders |19| . 
It is common to describe the quality of a road surface by a power law fit to the 
power spectrum of the height profile [2"Ull2 11122) . This method is not appropriate 
especially for wavelengths below 0.3 m and for non-Gaussian height distribu- 
tions. Furthermore, the power spectrum characterizes only the r-dependence 
of one moment of the two-point correlations . While some improvements to 
this method have been proposed [21] the characterization of road surfaces still 
remains incomplete |25H2fi| . For the improvement of road surface characteriza- 
tion multifractality and multiaffinity seem not to be appropriate tools because 
scaling is no constant feature of road surfaces. 

In the remainder of this letter we first introduce our method, the determi- 
nation of a Fokker-Planck equation for the evolution of conditional probability 
density functions (pdf) directly from experimental data. Next, we present a 
typical data set, show evidence of its Markov properties, and estimate the coef- 
ficients of the Fokker-Planck equation. At last, we evaluate the precision of the 
estimated coefficients by numerical reconstruction of conditional and uncondi- 
tional pdfs. 
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2 Method 



It is one common procedure to characterize the complexity of a rough surface 
by the statistics of the height increment |27| 

h r (x) := h(x + r/2) - h{x - r/2) (1) 

depending on the length scale r, as marked in fig. 2] Other scale dependent 
roughness measures can, for example, be found in |15I3| . Here we use the height 
increment h r because it is also directly linked to vehicle vibrations induced by 
the road surface if r is the wheelbase. Another argument for the use of h r 
is that its moments are connected with spatial correlation functions, but it 
should be pointed out that our method presented in the following could be 
easily generalized to any scale dependent measure, e.g. the above-mentioned 
w r (x). As a new ansatz, h r will be regarded here as a stochastic variable in 
r. Without loss of generality we consider the process as being directed from 
larger to smaller scales. Our interest is the investigation how surface roughness 
is linked between different length scales. 

Complete information about the stochastic process would be available by the 
knowledge of all possible n-point, or more precisely n-scale, probability density 
functions (pdf) p(h\, ri; h 2 , r 2 ; . . . ; h n , r n ) describing the probability of finding 
simultaneously the increments hi on the scale r±, hi on the scale r 2 , and so 
forth up to h n on the scale r n . Here we use the notation hi(x) — h ri (x), see Q. 
Without loss of generality we take r\ < r 2 < . . . < r n . As a first question one 
has to ask for a suitable simplification. In any case the n-scale joint pdf can be 
expressed by multiconditional pdf 

p(hi,r x ; . . .;h n ,r n ) = p(hi,ri\h 2 ,r 2 ; ■ ■ .;h n ,r n ) ■ p(h 2 ,r 2 \h 3 ,r 3 ; . . .;h n ,r n ) 

■ ... 'p(hn-i,r n -i\hn,r n ) -p{h n ,r n ) (2) 

where p(hi,ri\hj,rj) denotes a conditional probability, which is defined as the 
probability of finding the increment hi on the scale r» under the condition that 
simultaneously, i.e. at the same location x, on a larger scale rj the value hj was 
found. An important simplification arises if 

p(hi,n\hi+i,ri+v, ...;hn,r n )= p{h l ,r i \h l+ i,r l+ i) . (3) 

This property is the defining feature of a Markov process evolving from r,-+i to 
rj. Thus for a Markov process the n-scale joint pdf factorize into n conditional 
pdf 

p(hi,ri;...;h n ,r n ) =p(h 1 ,r 1 \h 2 ,r 2 ) ■ . . .-p(h 

)■ (4) 

This Markov property implies that the r-dependence of h r can be regarded as 
a stochastic process evolving in r, driven by deterministic and random forces. 
If additionally the included noise is Gaussian distributed, the process can be 
described by a Fokker-Planck equation j^H] . For our height profiles it takes the 
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form 



-r—p(h r ,r\h ,ro) = l-—D^(h r ,r) + -^D ( - 2 \h r ,r)^ p(h r ,r\h ,r ) . 

(5) 

The Fokker-Planck equation then describes the evolution of the conditional 
probability density function from larger to smaller length scales and thus also 
the complete n-scale statistics. The minus sign on the left side of eq. © ex- 
presses this direction of the process, furthermore the factor r corresponds to 
a logarithmic variable p = lnr which leads to simplified results in the case of 
scaling behaviour |29) . 

The term D^(h r ,r) is commonly denoted as drift term, describing the 
deterministic part of the process, and D^{h. r ,r) as diffusion term being the 
variance of a Gaussian, 5-correlated noise. Focus of our analysis is to show 
evidence of the above mentioned Markov property and to derive the drift and 
diffusion coefficients D^> and in eq. JSJ from experimental data. 



3 Data 




Figure 1: Cut-out from a height pro- 
file of an irregular cobblestone road. 
Additionally the construction of the 
height increment h r (xo) = h(xo + 
r/2) — h(xo — r/2) is illustrated. 
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Figure 2: Power spectral density of the 
data set. The arrow indicates begin- 
ning and direction of the range where 
drift and diffusion coefficient could be 
obtained and verified. 



Height profiles were measured from numerous road and cycle track surfaces 
typical for West German bicycle traffic. Road sections were selected in such a 
way that stationarity is given. Profile length is typically 20 m or 19 000 samples, 
respectively. Longitudinal resolution was 1.04 mm and vertical error smaller 
than 0.5 mm. As an example we present results from a data set of an irregular 
cobblestone road consisting of ten profiles with a total of about 190 000 sam- 
ples jHOl- Figure^Jshows a short section of the data. In fig- Elthe power spectral 
density of the height profiles is plotted against the wavenumber. Scaling be- 
haviour is not found at the beginning of the analysed region of length scales, 



while for smaller scales (20 < fc/m -1 < 100) it appears to be present. In the 
following, the height increments are normalized by := lim r ^ 00 (/i^) 1 ' 2 , with 
(Too = 6.3 mm for the given data set. 



4 Markov Properties 

For a Markov process the defining feature is that the n-point conditional proba- 
bility distributions are equal to the single conditional probabilities, according to 
eq. ©■ With the given amount of data points the verification of this condition 
is only possible for n — 3 and for n < r 2 < r?, < 300 mm. 
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Figure 3: a) Contour plot of single and double conditional probabilities 
p{hi,r 1 \h2,r 2 ) (dashed lines) and p(hi, r\\h,2, r%\ ha = 0, r^) (solid lines) for 
r\ = 126; r 2 = 188; r-s = 251 mm. Contour levels differ by a factor of 10, with 
an additional level at p = 0.3. b), c) Two one-dimensional cuts at h 2 ~ ±<7oo 
with p(h\, ri \h 2 , r 2 ) as dashed lines and p(hi, ri\h 2 , r 2 ; /13, r^j as circles. 



Figure |3 shows both sides of eq. © with n = 3 in a contour plot and two 
cuts at h 2 w itToo. The value of /13 was chosen to be /13 = 0. We take this rather 
good correspondence as a strong hint for a Markov process. Markov properties 
were found for scale distances from about 17 mm upwards. Note that the main 
axis of the distribution is tilted, indicating that p(h\,ri\h 2 ,r 2 ) ^ p(hi,ri) and 
thus height increments on different scales are not independent. 



5 Estimation of Drift and Diffusion Coefficients 

In order to obtain the drift (D^) and diffusion coefficient (D^) for eq. (J5J we 
proceed in a well defined way like it was already expressed by Kolmogorov |3J , 
see also |28II16| . First, the conditional moments M^ k '(h r , r, A r ) for finite step 
sizes Ar are directly estimated from the data via moments of the conditional 
probabilities 

M^(h r ,r,Ar) = -f— f + (h - h r f p(h,r - Ar\h r ,r) dh . (6) 
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Second, the coefficients (h r , r) are obtained from the limit o{M^(h r ,r, A r ) 
when Ar approaches zero 

D {k) (h r ,r)= limM(''(/i r ,r,Ar). (7) 




Figure 4: Estimated coefficients of the Fokker-Planck equation for r = 188 mm. 
Parameterizations are shown as solid and broken lines, a) Drift coefficient 
flWfV.r), b) Diffusion coefficient D^{h ri r). 

Figure 01 shows estimations of the drift coefficient and the diffusion 

coefficient for r = 188 mm. The error bars are estimated from the errors 
of M« and M< 2 ) via the number of events contributing to each value. The 
limit Ar — > was performed in both cases by a linear fit to M^ k \h r , r, Ar) 
in the range 17 mm < Ar < 29 mm. Both coefficients were parameterized 
as piecewise linear functions where the behaviour within \h r \ < 2.baoo could be 
derived directly from the above estimations. Outside this range increasing errors 
make a precise estimate more difficult. Parameterizations were chosen here 
to additionally obtain good results in the verification procedures (see below). 
Figure0]shows that the resulting parameterizations are in good agreement with 
the estimations. It is easy to verify that with linear and constant the 
Fokker-Planck equation J5J describes a Gaussian process, while with a parabolic 
Z)( 2 ) the distributions become non-Gaussian, also called intermittent or heavy 
tailed. 

In our case, here, D^(h r ,r) is characterized by the slope —j(r) within 
\h r \ < 2.5(700, — 4y(r) for h r < — 2.baoo, and —2^{r) for h r > 2.5(700 (compare 
fig. |3J|. The dependence of 7 on r is nontrivial with the value ranging from 
0.82 for r = 83mm to 1.9 for r = 188mm. D^(h r ,r) was found to have a 
value P(r) independent of h r within \h r \ < 2.5(7oo- For h r < — 2.5(7oo was 
parameterized to be linear with slope — 4/3(r), for h r > 2.5aoo with slope 3/3(r). 
The dependence of (3 on r can be approximated by (3{r) = 0.0117 r/mm. 



6 Verification of Coefficients 

Next, we want to evaluate the precision of our result. Therefore we return to 
eq. (JSJ. Knowing and it should be possible to calculate the pdf of h r 
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with the corresponding Fokker-Planck equation. Equation (jSJ) can be integrated 
over ho and is then valid also for the unconditional pdf. Now the empirical pdf 
at ro = 188 mm is parameterized (see fig. El and used as initial condition for 
a numerical solution of the integrated form of eq. ©. For several values of r 
the reconstructed pdf is compared to the respective empirical pdf, as shown in 
fig. |3J Please note that the interchange of steeper and flatter regions in the 
reconstructed pdf is achieved by the piecewise linear parameterization of D^ 
and DM. 




-4-3-2-1 1 2 3 4 

Figure 5: Numerical solution of the in- 
tegrated form of Fokker-Planck equa- 
tion (JjJ compared to empirical pdf 
(symbols) at different scales r. Solid 
line: empirical pdf parameterized at 
r = 188 mm, dashed lines: recon- 
structed pdf. Scales r are 188, 158, 
112, 79, 46 mm from top to bottom. 
Pdf are shifted in vertical direction for 
clarity of presentation. 




-3-1.5 1.5 3 -3-1.5 1.5 3 

h, KJ h, [oj 

Figure 6: Direct numerical solution of 
Fokker-Planck equation (J^J compared 
to the empirical pdf at scales ro = 188, 
7"i = 131 mm. a) Contour plot of em- 
pirical (solid lines) and reconstructed 
pdf (dashed lines). Contour levels are 
as in fig. 01 b,c) Cuts at ho ~ ic^. 
Empirical pdf are plotted as symbols. 



A second verification is the reconstruction of conditional pdf by direct nu- 
merical solution of the Fokker-Planck equation (J5J. An example for the scales 
ro = 188 mm and r% — 131mm is shown in fig. El Reconstructing the condi- 
tional pdf this way is much more sensitive to deviations in D^ and D^ . This 
becomes evident by the fact that the conditional pdf (and not the unconditional 
pdf of fig. [SJ) determine D^ and D^ (see © and J7J). Here also the differ- 
ence to the multiscaling analysis becomes clear, which analyses higher moments 
(hi.) = J h% ■ p(h r )dh r of h r , and does not depend on conditional pdf. It is 
easy to show that there are many different stochastic processes which lead to 
the same single scale pdf p(h r ). 
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7 Conclusions 



The height increment h r of surface height profiles as a stochastic variable in 
r can be correctly described by a Fokker-Planck equation with drift and diffu- 
sion coefficients derived directly from measured data. The results of the pre- 
sented example support the hypothesis that the noise term in the evolution 
of the stochastic variable h r in r is sufficiently well described by a Gaussian, 
5-correlated random process. 

As the Fokker-Planck equation describes the evolution oip(h r , r\ho,ro) and 
p(h r ,r) with r, it covers also the behaviour of the moments (h™) including any 
possible scaling behaviour. From the integrated form of eq. (JSJ) an equation for 
the moments can be obtained by multiplying with h™ and integrating over h r . 
For £)W being purely linear in h r and purely quadratic, multifractal scaling 
is obtained. We note again that, compared to scaling features, the knowledge of 
the Fokker-Planck equation provides more information on the complexity of the 
surface roughness in the sense of multi-scale joint probability density functions, 
eq. (j2J, which correspond to multipoint statistics. While to this end we do 
not seem to find universal laws concerning rough structures, we do achieve a 
comprehensive characterization of a specific surface, showing the strength and 
generality of this method. 

At last we want to point out that the Fokker-Planck equation corresponds 
to an equivalent Langevin equation |28j . The use of this Langevin equation in 
the scale variable should open the possibility to directly simulate surface profiles 
with given stochastic properties for different applications. 

Financial support by the Volkswagen Foundation is kindly acknowledged. 
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